﻿<!DOCTYPE html>
<html>
  <head>
    <title>Geodesic calculations for an ellipsoid done right</title>
    <meta name="description"
	  content="Geodesic calculations for an ellipsoid done right" />
    <meta name="author" content="Charles F. F. Karney">
    <meta name="keywords"
	  content="geodesics,
		   geodesic distance,
		   geographic distance,
		   shortest path,
		   direct geodesic problem,
		   inverse geodesic problem,
		   distance and azimuth,
		   distance and heading,
		   range and bearing,
		   geographic area,
		   geodesic polygon,
		   spheroidal triangle,
		   latitude and longitude,
		   online calculator,
		   WGS84 ellipsoid,
		   GeographicLib" />
    <script type="text/javascript"
	    src="http://geographiclib.sf.net/scripts/geographiclib.js">
    </script>
    <script type="text/javascript">
<!--
var geod = GeographicLib.Geodesic.WGS84;
var dms = GeographicLib.DMS;
function formatpoint(lat, lon, azi, dmsformat, prec) {
  prec += 5;
  if (dmsformat) {
    var trail = prec < 2 ? dms.DEGREE :
      (prec < 4 ? dms.MINUTE : dms.SECOND);
    prec = prec < 2 ? prec : (prec < 4 ? prec - 2 : prec - 4);
    return (dms.Encode(lat, trail, prec, dms.LATITUDE) + " " +
	    dms.Encode(lon, trail, prec, dms.LONGITUDE) + " " +
	    dms.Encode(azi, trail, prec, dms.AZIMUTH));
  } else {
    return (lat.toFixed(prec) + " " +
	    lon.toFixed(prec) + " " +
	    azi .toFixed(prec));
  }
}

function GeodesicInverse(input, dmsformat, prec) {
  var result = {};
  try {
    // Input is a blank-delimited line: lat1 lon1 lat2 lon2
    var t = new String(input);
    t = t.replace(/^\s+/,"").replace(/\s+$/,"").split(/[\s,]+/,6);
    if (t.length != 4)
      throw new Error("Need 4 input items");
    var p1 = GeographicLib.DMS.DecodeLatLon(t[0], t[1]);
    var p2 = GeographicLib.DMS.DecodeLatLon(t[2], t[3]);
    t = geod.Inverse(p1.lat, p1.lon, p2.lat, p2.lon);
    result.status = "OK";
    result.p1 = formatpoint(t.lat1, t.lon1, t.azi1, dmsformat, prec);
    result.p2 = formatpoint(t.lat2, t.lon2, t.azi2, dmsformat, prec);
    result.s12 = t.s12.toFixed(prec);
  }
  catch (e) {
    result.status = "ERROR: " + e.message;
    result.p1 = "";
    result.p2 = "";
    result.s12 = "";
  }
  return result;
}

function GeodesicDirect(input, dmsformat, prec) {
  var result = {};
  try {
    // Input is a blank-delimited line: lat1 lon1 azi1 s12
    var t = new String(input);
    t = t.replace(/^\s+/,"").replace(/\s+$/,"").split(/[\s,]+/,6);
    if (t.length != 4)
      throw new Error("Need 4 input items");
    var p1 = GeographicLib.DMS.DecodeLatLon(t[0], t[1]);
    var azi1 = GeographicLib.DMS.DecodeAzimuth(t[2]);
    var s12 = parseFloat(t[3]);
    t = geod.Direct(p1.lat, p1.lon, azi1, s12);
    result.status = "OK";
    result.p1 = formatpoint(t.lat1, t.lon1, t.azi1, dmsformat, prec);
    result.p2 = formatpoint(t.lat2, t.lon2, t.azi2, dmsformat, prec);
    result.s12 = t.s12.toFixed(prec);
  }
  catch (e) {
    result.status = "ERROR: " + e.message;
    result.p1 = "";
    result.p2 = "";
    result.s12 = "";
  }
  return result;
}

function GeodesicInversePath(input, dmsformat, prec) {
  var result = {};
  try {
    // Input is a blank-delimited line: lat1 lon1 lat2 lon2 ds12 maxnum
    var t = new String(input);
    t = t.replace(/^\s+/,"").replace(/\s+$/,"").split(/[\s,]+/,8);
    if (t.length != 6)
      throw new Error("Need 6 input items");
    var p1 = GeographicLib.DMS.DecodeLatLon(t[0], t[1]);
    var p2 = GeographicLib.DMS.DecodeLatLon(t[2], t[3]);
    var ds12 = parseFloat(t[4]);
    var maxnum = parseInt(t[5]);
    var t = new String(input);
    t = t.split(/[\s,]+/,8);
    if (t[0] == "") t.shift();
    t = geod.InversePath(p1.lat, p1.lon, p2.lat, p2.lon, ds12, maxnum);
    result.status = "OK";
    result.points = ""
    for (var i = 0; i < t.length; ++i)
      result.points +=
      formatpoint(t[i].lat, t[i].lon, t[i].azi, dmsformat, prec) + "\n";
  }
  catch (e) {
    result.status = "ERROR: " + e.message;
    result.points = "";
  }
  return result;
}

function GeodesicArea(input, polyline) {
  var result = {};
  try {
    // Input is a newline-delimited points;
    // each point is blank-delimited: lat lon
    var t = new String(input);
    t = t.split(/[\r\n]/);
    if (t[0] == "") t.shift();
    if (t[t.length-1] == "") t.pop();
    for (var i = 0; i < t.length; ++i) {
      var pos = t[i].replace(/^\s+/,"").replace(/\s+$/,"").split(/[\s,]+/,4);
      if (pos.length != 2)
	throw new Error("Need 2 items on each line");
      t[i] = dms.DecodeLatLon(pos[0], pos[1]);
    }
    t = geod.Area(t, polyline);
    result.status = "OK";
    result.area = t.number + " " +
      t.perimeter.toFixed(8);
    if (!polyline)
      result.area += " " + t.area.toFixed(2);
  }
  catch (e) {
    result.status = "ERROR: " + e.message;
    result.area = "";
  }
  return result;
}
//-->
    </script>
  </head>
  <body>
    <div>
      <h2>Geodesic calculations for an ellipsoid done right</h2>
    </div>
    <p>
      This page illustrates the geodesic routines available in
      <a href="http://geographiclib.sourceforge.net">GeographicLib</a>.
      The C++ code has been converted to JavaScript so the calculations
      are carried out on the client.  The algorithms are considerably
      more accurate than Vincenty's method, and offer more functionality
      (an inverse method which never fails to converge, differential
      properties of the geodesic, and the area under a geodesic).  The
      algorithms are derived in
      <blockquote>
	Charles F. F. Karney,<br>
	<a href="https://dx.doi.org/10.1007/s00190-012-0578-z">
	  <i>Algorithms for geodesics</i></a>,<br>
	J. Geodesy <b>87</b>(1), 43&ndash;55 (Jan. 2013);<br>
	DOI:
	<a href="https://dx.doi.org/10.1007/s00190-012-0578-z">
	  10.1007/s00190-012-0578-z</a>
	(<a href="https://dx.doi.org/10.1007/s00190-012-0578-z">pdf</a>);
	addenda: <a href="http://geographiclib.sf.net/geod-addenda.html">
	  geod-addenda.html</a>.
      </blockquote>
      This page just provides a basic interface.  Enter latitudes,
      longitudes, and azimuths as degrees and distances as meters using
      spaces or commas as separators.  (Angles may be entered as decimal
      degrees or as degrees, minutes, and seconds, e.g. -20.51125,
      20&deg;30&prime;40.5&Prime;S, S20d30'40.5&quot;, or
      -20:30:40.5.)  The results are accurate to about
      15&nbsp;nanometers (or 0.1&nbsp;m<sup>2</sup> per vertex for
      areas).  A slicker page where the geodesics are incorporated into
      Google Maps is given <a href="geod-google.html">here</a>.
    <p>
      Jump to:
      <ul>
	<li><a href="#inverse">Inverse problem</a></li>
	<li><a href="#direct">Direct problem</a></li>
	<li><a href="#path">Geodesic path</a></li>
	<li><a href="#area">Polygon area</a></li>
	<li><a href="geod-google.html">Geodesic lines, circles, and
	envelopes in Google Maps</a></li>
      </ul>
    </p>
    <hr>
    <form name="inverse" >
      <h3><a class="anchor" id="inverse">Inverse problem</h3>
      <p>
	Find the shortest path between two points on the earth.  The
	path is characterized by its length <i>s12</i> and its azimuth
	at the two ends <i>azi1</i> and <i>azi2</i>.  The sample
	calculation finds the shortest path between Wellington, New
	Zealand, and Salamanca, Spain.  (For this example, the
	<a href="http://www.ngs.noaa.gov/">NGS</a>
	<a href="http://www.ngs.noaa.gov/cgi-bin/Inv_Fwd/inverse2.prl">
        inverse geodesic calculator</a> returns a result which is 1.2 km
	too long with an azimuth which is off by 3 degrees.)  To perform
	the calculation, press the &ldquo;COMPUTE&rdquo; button.
      </p>
      <p>Enter <i>&ldquo;lat1 lon1 lat2 lon2&rdquo;</i>:</p>
      <p>input:
        <input name="input" size=72 value="-41.32 174.81 40.96 -5.50" />
      </p>
      <p>
	Output format:&nbsp;&nbsp;<label for="ig">
	  <input type="radio" value="g" name="format" id="ig" checked>
	  decimal degrees
	</label>&nbsp;
	<label for="id">
	  <input type="radio" value="d" name="format" id="id">
	  degrees minutes seconds
	</label><br>
	Output precision:&nbsp;&nbsp;<select name="prec" size=1>
	  <option value='0'> 1m 0.00001d 0.1"</option>
	  <option value='1'> 100mm 0.01"</option>
	  <option value='2'> 10mm 0.001"</option>
	  <option value='3' selected> 1mm 0.0001"</option>
	  <option value='4'> 100um 0.00001"</option>
	  <option value='5'> 10um 0.000001"</option>
	  <option value='6'> 1um 0.0000001"</option>
	  <option value='7'> 100nm 0.00000001"</option>
	  <option value='8'> 10nm 0.000000001"</option>
	  <option value='9'> 1nm 0.0000000001"</option>
	</select>
      </p>
      <p>
        <input type="button" value="COMPUTE"
	       onclick="var t = GeodesicInverse(document.inverse.input.value,
			        document.inverse.format[1].checked,
				document.inverse.prec.selectedIndex);
			document.inverse.status.value = t.status;
			document.inverse.p1.value = t.p1;
			document.inverse.p2.value = t.p2;
			document.inverse.s12.value = t.s12;" />
      </p>
      <p>
	status:
        <input name="status" size=50 readonly />
      </p>
      <p>
	lat1 lon1 <font color='blue'>azi1</font>:
        <input name="p1" size=75 readonly />
      </p>
      <p>
	lat2 lon2 <font color='blue'>azi2</font>:
        <input name="p2" size=75 readonly />
      </p>
      <p>
	<font color='blue'>s12</font>:
        <input name="s12" size=25 readonly />
      </p>
    </form>
    <hr>
    <form name="direct">
      <h3><a class="anchor" id="direct">Direct problem</h3>
      <p>
	Find the destination traveling a given distance along a geodesic
	with a given azimuth at the starting point.  The destination is
	characterized by its position <i>lat2, lon2</i> and its azimuth
	at the destination <i>azi2</i>.  The sample calculation shows
	the result of travelling 10000 km NE from JFK airport.  To perform
	the calculation, press the &ldquo;COMPUTE&rdquo; button.
      </p>
      <p>Enter <i>&ldquo;lat1 lon1 azi1 s12&rdquo;</i>:</p>
      <p>input:
        <input name="input" size=72 value="40.6 -73.8 45 10000e3" />
      </p>
      <p>
	Output format:&nbsp;&nbsp;<label for="dg">
	  <input type="radio" value="g" name="format" id="dg" checked>
	  decimal degrees
	</label>&nbsp;
	<label for="dd">
	  <input type="radio" value="d" name="format" id="dd">
	  degrees minutes seconds
	</label><br>
	Output precision:&nbsp;&nbsp;<select name="prec" size=1>
	  <option value='0'> 1m 0.00001d 0.1"</option>
	  <option value='1'> 100mm 0.01"</option>
	  <option value='2'> 10mm 0.001"</option>
	  <option value='3' selected> 1mm 0.0001"</option>
	  <option value='4'> 100um 0.00001"</option>
	  <option value='5'> 10um 0.000001"</option>
	  <option value='6'> 1um 0.0000001"</option>
	  <option value='7'> 100nm 0.00000001"</option>
	  <option value='8'> 10nm 0.000000001"</option>
	  <option value='9'> 1nm 0.0000000001"</option>
	</select>
      </p>
      <p>
        <input type="button" value="COMPUTE"
	       onclick="var t = GeodesicDirect(document.direct.input.value,
			        document.direct.format[1].checked,
				document.direct.prec.selectedIndex);
			document.direct.status.value = t.status;
			document.direct.p1.value = t.p1;
			document.direct.p2.value = t.p2;
			document.direct.s12.value = t.s12;" />
      </p>
      <p>
	status:
        <input name="status" size=50 readonly />
      </p>
      <p>
	lat1 lon1 azi1:
        <input name="p1" size=75 readonly />
      </p>
      <p>
	<font color='blue'>lat2 lon2 azi2</font>:
        <input name="p2" size=75 readonly />
      </p>
      <p>
	s12:
        <input name="s12" size=25 readonly />
      </p>
    </form>
    <hr>
    <form name="path">
      <h3><a class="anchor" id="path">Geodesic path</h3>
      <p>
	Find intermediate points along a geodesic.  In addition to
	specifying the endpoints, give <i>ds12</i>, the maximum distance
	between the intermediate points and <i>maxk</i>, the maximum
	number of intervals the geodesic is broken into.  The output
	gives a sequence of positions <i>lat, lon</i> together with the
	corresponding azimuths <i>azi</i>.  The sample shows the path
	from JFK to Singapore's Changi Airport at about 1000 km
	intervals.  (In this example, the path taken by Google Earth
	deviates from the shortest path by about 2.9 km.)  To perform
	the calculation, press the &ldquo;COMPUTE&rdquo; button.
      </p>
      <p>Enter <i>&ldquo;lat1 lon1 lat2 lon2 ds12 maxk&rdquo;</i>:</p>
      <p>input:
        <input name="input" size=72 value="40.6 -73.8 1.4 104 1000e3 20" />
      </p>
      <p>
	Output format:&nbsp;&nbsp;<label for="pg">
	  <input type="radio" value="g" name="format" id="pg" checked>
	  decimal degrees
	</label>&nbsp;
	<label for="pd">
	  <input type="radio" value="d" name="format" id="pd">
	  degrees minutes seconds
	</label><br>
	Output precision:&nbsp;&nbsp;<select name="prec" size=1>
	  <option value='0' selected> 1m 0.00001d 0.1"</option>
	  <option value='1'> 100mm 0.01"</option>
	  <option value='2'> 10mm 0.001"</option>
	  <option value='3'> 1mm 0.0001"</option>
	  <option value='4'> 100um 0.00001"</option>
	  <option value='5'> 10um 0.000001"</option>
	  <option value='6'> 1um 0.0000001"</option>
	  <option value='7'> 100nm 0.00000001"</option>
	  <option value='8'> 10nm 0.000000001"</option>
	  <option value='9'> 1nm 0.0000000001"</option>
	</select>
      </p>
      <p>
        <input type="button" value="COMPUTE"
	       onclick="var t = GeodesicInversePath(document.path.input.value,
			        document.path.format[1].checked,
				document.path.prec.selectedIndex);
			document.path.status.value = t.status;
			document.path.points.value = t.points;" />
      </p>
      <p>
	status:
        <input name="status" size=50 readonly />
      </p>
      <p>
	points (lat lon azi):<br>
        <textarea name="points" cols=70 rows=21 readonly></textarea>
      </p>
    </form>
    <hr>
    <form name="area">
      <h3><a class="anchor" id="area">Polygon area</h3>
      <p>
	Find the perimeter and area of a polygon whose sides are
	geodesics.  The polygon must be simple (i.e., must not intersect
	itself).  (There's no need to ensure that the polygon is
	closed.)  Counter-clockwise traversal of the polygon results in
	a positive area.  The polygon can encircle one or both poles.
	The sample gives the approximate perimeter (in m) and area (in
	m<sup>2</sup>) of Antarctica.  (For this example, Google Earth
	Pro returns an area which is 30 times too large!  However this
	is a little unfair, since Google Earth has no concept of
	polygons which encircle a pole.)  If the <i>polyline</i> option
	is selected then just the length of the line joining the points
	is returned.  To perform the calculation, press the
	&ldquo;COMPUTE&rdquo; button.
      </p>
      <p>Enter points, one per line, as <i>&ldquo;lat lon&rdquo;</i>:</p>
      <p>points (lat lon):<br>
        <textarea name="input" cols=36 rows=13>-63.1  -58
-72.9  -74
-71.9 -102
-74.9 -102
-74.3 -131
-77.5 -163
-77.4  163
-71.7  172
-65.9  140
-65.7  113
-66.6   88
-66.9   59
-69.8   25
-70.0   -4
-71.0  -14
-77.3  -33
-77.9  -46
-74.7  -61
</textarea>
      </p>
      <p>
	Treat points as:&nbsp;&nbsp;<label for="lg">
	  <input type="radio" value="p" name="polyline" id="lg" checked>
	  polygon
	</label>&nbsp;
	<label for="ll">
	  <input type="radio" value="l" name="polyline" id="ll">
	  polyline
	</label>
      </p>
      <p>
        <input type="button" value="COMPUTE"
	       onclick="var t = GeodesicArea(document.area.input.value,
			document.area.polyline[1].checked);
			document.area.status.value = t.status;
			document.area.area.value = t.area;" />
      </p>
      <p>
	status:
        <input name="status" size=50 readonly />
      </p>
      <p>
	number perimeter area:
        <input name="area" size=55 readonly />
      </p>
    </form>
    <hr>
    <address>Charles Karney
      <a href="mailto:charles@karney.com">&lt;charles@karney.com&gt;</a>
      (2011-08-04)</address>
    <br>
    <a href="http://geographiclib.sourceforge.net">Geographiclib Sourceforge</a>
  </body>
</html>
